#include <fstream>
#include <iostream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
#include "TLine.h"
#include "TSystem.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TArc.h"
#include "TEllipse.h"
#include "TClonesArray.h"

using namespace std;

void Ar(){

  int NO = 10;

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

  TFile *fread=new TFile("/scratch/shinichi/OutputFiles/test.root");

  if(!fread){ cout<<"No Such File"<<endl; return;}

  TTree *tree=(TTree*)fread->Get("tree");
  
  TClonesArray *JPMTa=new TClonesArray("PMT");
  TClonesArray *MPMT1a=new TClonesArray("PMT");
  TClonesArray *MPMT2a=new TClonesArray("PMT");
  TClonesArray *MPMT3a=new TClonesArray("PMT");
  TClonesArray *MPMT4a=new TClonesArray("PMT");
  TClonesArray *Particlea=new TClonesArray("Particle");
  TClonesArray *Neutrona=new TClonesArray("Particle");
  TClonesArray *Argona=new TClonesArray("Particle");
  TClonesArray *Gammaa=new TClonesArray("Particle");
  int nPhoton,nElectron;
  int nJhit,nM1hit,nM2hit,nM3hit,nM4hit;

  tree->GetBranch("JPMTb")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT1b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT2b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT3b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT4b")->SetAutoDelete(kFALSE);
  tree->GetBranch("Particleb")->SetAutoDelete(kFALSE);
  tree->SetBranchAddress("JPMTb",&JPMTa);
  tree->SetBranchAddress("MPMT1b",&MPMT1a);
  tree->SetBranchAddress("MPMT2b",&MPMT2a);
  tree->SetBranchAddress("MPMT3b",&MPMT3a);
  tree->SetBranchAddress("MPMT4b",&MPMT4a);
  tree->SetBranchAddress("Particleb",&Particlea);
  tree->SetBranchAddress("Neutronb",&Neutrona);
  tree->SetBranchAddress("Argonb",&Argona);
  tree->SetBranchAddress("Gammab",&Gammaa);
  tree->SetBranchAddress("nPhoton",&nPhoton);
  tree->SetBranchAddress("nElectron",&nElectron);
  tree->SetBranchAddress("nJPMThits",&nJhit);
  tree->SetBranchAddress("nMPMT1hits",&nM1hit);
  tree->SetBranchAddress("nMPMT2hits",&nM2hit);
  tree->SetBranchAddress("nMPMT3hits",&nM3hit);
  tree->SetBranchAddress("nMPMT4hits",&nM4hit);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;

  TGraph *ArgonXY[32];
  TGraph *ArgonXZ[32];
  TGraph *GammaXY[32];
  TGraph *GammaXZ[32];
  for(int i=0;i<32;i++){
    ArgonXY[i]=new TGraph();;
    ArgonXZ[i]=new TGraph();;
    GammaXY[i]=new TGraph();
    GammaXZ[i]=new TGraph();
  }

  TH1F *hNPhoton[8];
  TH1F *hNElectron[8];
  for(int i=0;i<8;i++){
    hNPhoton[i]=new TH1F(Form("hPhoton%d",i),Form("hPhoton%d",i),1000,0,100000);
    hNElectron[i]=new TH1F(Form("hElectron%d",i),Form("hElectron%d",i),500,0,50000);
  }
  TH2F *hJPMTxy=new TH2F("hJPMTxy","hJPMTxy",150,-150,150,150,-150,150);
  TH2F *hMPMT1xy=new TH2F("hMPMT1xy","hMPMT1xy",150,-150,150,150,-150,150);
  TH2F *hMPMT2xy=new TH2F("hMPMT2xy","hMPMT2xy",150,-150,150,150,-150,150);
  TH2F *hMPMT3xy=new TH2F("hMPMT3xy","hMPMT3xy",150,-150,150,150,-150,150);
  TH2F *hMPMT4xy=new TH2F("hMPMT4xy","hMPMT4xy",150,-150,150,150,-150,150);

  TH2F   *hPMTeff=new TH2F("hPMTeff","hPMTeff",30,-150,150,35,-650,-300);
  TH2F  *hJPMTeff=new TH2F("hJPMTeff","hJPMTeff",30,-150,150,35,-650,-300);
  TH2F *hM1PMTeff=new TH2F("hM1PMTeff","hM1PMTeff",30,-150,150,35,-650,-300);
  TH2F *hM2PMTeff=new TH2F("hM2PMTeff","hM2PMTeff",30,-150,150,35,-650,-300);
  TH2F *hM3PMTeff=new TH2F("hM3PMTeff","hM3PMTeff",30,-150,150,35,-650,-300);
  TH2F *hM4PMTeff=new TH2F("hM4PMTeff","hM4PMTeff",30,-150,150,35,-650,-300);

  int NARGON,NGAMMA;

  for(int iev=0;iev<nev;iev++){
    JPMTa->Clear();
    MPMT1a->Clear();
    MPMT2a->Clear();
    MPMT3a->Clear();
    MPMT4a->Clear();
    Particlea->Clear();
    Neutrona->Clear();
    Argona->Clear();
    Gammaa->Clear();

    tree->GetEntry(iev);

    //----- Argon
    int iargon = 0;
    int trackid = -1;
    int istep = 0;
    for(int i=0;i<Argona->GetEntriesFast();i++){
      Particle *argon = (Particle*)Argona->At(i);
      if(trackid!=argon->TrackID){
	trackid = argon->TrackID;
	istep = 0;
	iargon++;
	if( iev==NO ){
	  ArgonXY[iargon-1]->SetPoint(istep,argon->xpre,argon->ypre);
	  ArgonXZ[iargon-1]->SetPoint(istep,argon->xpre,argon->zpre);
	}
      }
      if( iev==NO ){
	ArgonXY[iargon-1]->SetPoint(istep+1,argon->xpos,argon->ypos);
	ArgonXZ[iargon-1]->SetPoint(istep+1,argon->xpos,argon->zpos);
      }
      istep++;
    }
    if( iev==NO ){
      NARGON = iargon;
      cout<<"nArgon = "<<NARGON<<endl;
    }

    //----- Gamma
    int igamma = 0;
    trackid=-1;
    istep = 0;
    for(int i=0;i<Gammaa->GetEntriesFast();i++){
      Particle *gamma = (Particle*)Gammaa->At(i);
      if(trackid!=gamma->TrackID){
	trackid = gamma->TrackID;
	istep = 0;
	igamma++;
	if( iev==NO ){
	  GammaXY[igamma-1]->SetPoint(istep,gamma->xpre,gamma->ypre);
	  GammaXZ[igamma-1]->SetPoint(istep,gamma->xpre,gamma->zpre);
	}
      }
      if( iev==NO ){
	//	cout<<istep<<" edep = "<<gamma->edep<<endl;
	GammaXY[igamma-1]->SetPoint(istep+1,gamma->xpos,gamma->ypos);
	GammaXZ[igamma-1]->SetPoint(istep+1,gamma->xpos,gamma->zpos);
      }
      istep++;
    }
    if( iev==NO ){
      NGAMMA = igamma;
      cout<<"nGamma = "<<NGAMMA<<endl;
    }

    //----- PMT efficiency
    int iz = iev / 21 % 18;
    int ir = iev % 21;
    float zpos = iz*10.-545;
    float rpos = ir*10.-120;
    hPMTeff->Fill( rpos,zpos,(nJhit+nM1hit+nM2hit+nM3hit+nM4hit)/(float)nPhoton );
     hJPMTeff->Fill( rpos,zpos, nJhit/(float)nPhoton );
    hM1PMTeff->Fill( rpos,zpos,nM1hit/(float)nPhoton );
    hM2PMTeff->Fill( rpos,zpos,nM2hit/(float)nPhoton );
    hM3PMTeff->Fill( rpos,zpos,nM3hit/(float)nPhoton );
    hM4PMTeff->Fill( rpos,zpos,nM4hit/(float)nPhoton );

    //----- PMT Hit
    if( iev==NO ){
      for(int iJhit=0;iJhit<nJhit;iJhit++){
	PMT *jpmt = (PMT*)JPMTa->At(iJhit);
	hJPMTxy->Fill( jpmt->xpos,jpmt->ypos );
      }
      for(int iM1hit=0;iM1hit<nM1hit;iM1hit++){
	PMT *mpmt1 = (PMT*)MPMT1a->At(iM1hit);
	hMPMT1xy->Fill( mpmt1->xpos,mpmt1->ypos );
      }
      for(int iM2hit=0;iM2hit<nM2hit;iM2hit++){
	PMT *mpmt2 = (PMT*)MPMT2a->At(iM2hit);
	hMPMT2xy->Fill( mpmt2->xpos,mpmt2->ypos );
      }
      for(int iM3hit=0;iM3hit<nM3hit;iM3hit++){
	PMT *mpmt3 = (PMT*)MPMT3a->At(iM3hit);
	hMPMT3xy->Fill( mpmt3->xpos,mpmt3->ypos );
      }
      for(int iM4hit=0;iM4hit<nM4hit;iM4hit++){
	PMT *mpmt4 = (PMT*)MPMT4a->At(iM4hit);
	hMPMT4xy->Fill( mpmt4->xpos,mpmt4->ypos );
      }
    }


  }

  //----- Draw Part
  //___ Event Display
  TCanvas *c1=new TCanvas("c1","c1",1200,600);
  c1->SetGrid();
  c1->Divide(2,1);
  c1->cd(1);
  TH1 *frame1=gPad->DrawFrame(-180,-180,180,180);
  frame1->SetXTitle("X (mm)");
  frame1->SetYTitle("Y (mm)");
  frame1->Draw();
  c1->cd(2);
  TH1 *frame2=gPad->DrawFrame(-180,-650,180,-300);
  frame2->SetXTitle("X (mm)");
  frame2->SetYTitle("Z (mm)");
  frame2->Draw();

  for(int igamma=0;igamma<NGAMMA;igamma++){
    c1->cd(1);
    GammaXY[igamma]->SetLineColor(4);
    GammaXY[igamma]->SetMarkerColor(4);
    GammaXY[igamma]->SetMarkerStyle(21);
    GammaXY[igamma]->SetMarkerSize(1);
    GammaXY[igamma]->Draw("pl");
    c1->cd(2);
    GammaXZ[igamma]->SetLineColor(4);
    GammaXZ[igamma]->SetMarkerColor(4);
    GammaXZ[igamma]->SetMarkerStyle(21);
    GammaXZ[igamma]->SetMarkerSize(1);
    GammaXZ[igamma]->Draw("pl");
  }
  for(int iargon=0;iargon<NARGON;iargon++){
    c1->cd(1);
    ArgonXY[iargon]->SetLineColor(3);
    ArgonXY[iargon]->SetMarkerColor(3);
    ArgonXY[iargon]->SetMarkerStyle(21);
    ArgonXY[iargon]->SetMarkerSize(1);
    ArgonXY[iargon]->Draw("pl");
    c1->cd(2);
    ArgonXZ[iargon]->SetLineColor(3);
    ArgonXZ[iargon]->SetMarkerColor(3);
    ArgonXZ[iargon]->SetMarkerStyle(21);
    ArgonXZ[iargon]->SetMarkerSize(1);
    ArgonXZ[iargon]->Draw("pl");
  }

  TArc *aInner=new TArc(0.,0.,133.,45.,315.);
  TLine *lInner=new TLine(133/sqrt(2),-133/sqrt(2),133/sqrt(2),133/sqrt(2));
  aInner->SetFillStyle(0);
  aInner->SetNoEdges(1);
  aInner->Draw();
  lInner->Draw();
  TLine *lFid1=new TLine(-133,-360,-133,-550);
  TLine *lFid2=new TLine(-133,-550,  94,-550);
  TLine *lFid3=new TLine(  94,-550,  94,-360);
  TLine *lFid4=new TLine(  94,-360,-133,-360);

  c1->cd(1);
  aInner->Draw();
  lInner->Draw();
  c1->cd(2);
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();

  //___ Photon Distribution
  /*
  TCanvas *c2=new TCanvas("c2","c2",800,600);
  hJPMTxy->Draw("colz");
  hMPMT1xy->Draw("samecolz");
  hMPMT2xy->Draw("samecolz");
  hMPMT3xy->Draw("samecolz");
  hMPMT4xy->Draw("samecolz");
  */

  //___ Photon eff
  gStyle->SetOptStat(000000);
  TCanvas *c3=new TCanvas("c3","c3",1200,900);
  c3->Divide(3,3);
  c3->cd(5);
  hJPMTeff->Draw("colz");
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();
  c3->cd(4);
  hM1PMTeff->Draw("colz");
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();
  c3->cd(8);
  hM2PMTeff->Draw("colz");
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();
  c3->cd(6);
  hM3PMTeff->Draw("colz");
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();
  c3->cd(2);
  hM4PMTeff->Draw("colz");
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();




}
